# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

from fractions import Fraction
from copy import copy

def divide_polynomials(fs, g):
    coefficients = [Polynomial([(Fraction(0), (0, 0))]) for x in fs]
    r = Polynomial()
    p = g

    while not zero_polynomial(p):
        i = 1
        divisionoccurred = False

        while i <= len(fs) and divisionoccurred != True:
            if monomial_divisor(p.leading_term(), fs[i-1].leading_term()):
                coefficients[i-1] = add_polynomials(
                    coefficients[i-1], 
                    Polynomial(
                        [divide_monomials(
                                p.leading_term(), fs[i-1].leading_term())]))

                p = add_polynomials(
                    p, multiply_polynomials(
                        Polynomial(
                            [divide_monomials(
                                    p.leading_term(), 
                                    fs[i-1].leading_term())]), 
                        fs[i-1]).negate())
                divisionoccurred = True
            else:
                i = i + 1
        if divisionoccurred != True:
            r = add_polynomials(r, Polynomial([p.leading_term()]))
            p = add_polynomials(p, Polynomial([p.leading_term()]).negate())

    return (coefficients, r)

class Polynomial:
    def __init__(self, f=[], order="lex"):
        # TODO(patrick): Accept other orders to instantiate a Polynomial.
        if order != "lex":
            raise Exception("Order not implemented!")
        self.order = order

        # Ensure that all elements of f are Monomials, and have access to the
        # Monomial methods. We could spare those elements that are already a
        # Monomial from this computation, however that is unneccessary
        # complexity for an application of this size.
        # The addition of the next few lines allows one to create a Polynomial
        # by hand, whilst keeping the Polynomial methods suitable for use by
        # other methods such as add_polynomials().
        if not all(map(lambda x: isinstance(x, Monomial), f)):
            raw_monomials = map(lambda x: x.__repr__(), f)
            if isinstance(raw_monomials, str):
                raw_monomials = eval(raw_monomials[0])
            f = map(lambda x: Monomial(eval(x)), raw_monomials)

        # Explicitly copying monomials from f to instanciated object
        self.monomials = [monomial for monomial in f]

        # Many other approaches are equivalent:
        #   self.monomials = f
        # or
        #   self.monomials = []
        #   for monomial in f:
        #     self.monomials.append(monomial)
        # I find the for loop and the list comprehension better because
        # it emphasizes the notion that f is a list of monomials,
        # rather than a black box object.
    def add_monomial(self, f):
        self.monomials.append(f)
    def leading_term(self):
        leading_term = [0]*max(map(
            lambda x: len(x),
            [monomial.degrees for monomial in self.monomials]))
        leading_coefficient = Fraction(0)
        for monomial in self.monomials:
            if max(monomial.degrees) == 0 and max(leading_term) == 0:
                leading_coefficient = monomial.coefficient
                continue
            if max(monomial.degrees) < max(leading_term):
                continue
            elif max(monomial.degrees) >= max(leading_term):
                skip = False
                if max(monomial.degrees) == max(leading_term):
                    for i in range(len(monomial.degrees)):
                        if monomial.degrees[i] > leading_term[i]:
                            Skip = False
                            break
                        skip = True

                if skip == True:
                    continue

                leading_term = [0 for x in leading_term]
                for j in range(len(monomial.degrees)):
                    leading_term[j] += monomial.degrees[j]
                    leading_coefficient = monomial.coefficient
                continue

        return Monomial((leading_coefficient, tuple(leading_term)))
    def __repr__(self):
        return str([monomial.__repr__() for monomial in self.monomials])
    def negate(self):
        for monomial in self.monomials:
            monomial.negate()

        return self
        
class Monomial:
    # TODO(patrick): Accept other orders to instantiate a Monomial.
    def __init__(self, f, order="lex"):
        #if order != "lex":
        #    raise Exception("Order not implemented!")
        self.order = order

        if not isinstance(f[0], Fraction):
            raise Exception("Coefficient must be a Fraction!")
        self.coefficient = f[0]

        if not isinstance(f[1], tuple):
            raise Exception("Degrees must be contained in a tuple!")
        if not all(map(lambda x: isinstance(x, int), f[1])):
            raise Exception("Degrees must be ints.")
        self.degrees = f[1]
    def __repr__(self):
        return str((self.coefficient, self.degrees))
    def negate(self):
        self.coefficient = -self.coefficient

def same_monomial_degree(f, g):
    if f.degrees == g.degrees:
        return True
    else:
        return False

def add_monomials(f, g):
    if f.degrees != g.degrees:
        raise Exception("Monomials are not of the same degree!")

    coefficient = f.coefficient + g.coefficient
    return Monomial((coefficient, f.degrees))

def multiply_monomials(f, g):
    check_ordering([f, g])

    coefficient = f.coefficient * g.coefficient

    degrees = [f.degrees, g.degrees]

    new_degrees = [0]*max(map(lambda x: len(x), degrees))

    for i in range(len(new_degrees)):
        for degree_tuple in degrees:
            try:
                new_degrees[i] += degree_tuple[i]
            except IndexError, e:
                continue

    return Monomial((coefficient, tuple(new_degrees)))

def monomial_divisor(f, g):
    for i in range(len(f.degrees)):
        if g.degrees[i] > f.degrees[i]:
            return False
    
    return True

def divide_monomials(f, g):
    if not monomial_divisor(f, g):
        return None

    coefficient = f.coefficient / g.coefficient

    degrees = [0]*len(f.degrees)
    for i in range(len(f.degrees)):
        degrees[i] = f.degrees[i] - g.degrees[i]

    return Monomial((coefficient, tuple(degrees)))

def add_polynomials(f, g):
    check_ordering([f, g])

    h = Polynomial()

    # Copy objects, so that we aren't modifying the original objects
    g_monomials = copy(g.monomials)
    f_monomials = copy(f.monomials)

    for monomial in f_monomials:
        # Remove coefficients from monomials
        degrees = map(lambda x: x.degrees, g_monomials)

        new_monomial = None
        for degree in degrees:
            # Monomial of matching degree found
            if monomial.degrees == degree:
                # Get original instances of monomials
                f_monomial = monomial
                g_monomial = g_monomials[degrees.index(degree)]

                # Remove matching monomial from list of monomials to enumerate
                g_monomials.remove(g_monomials[degrees.index(degree)])
                
                new_monomial = add_monomials(f_monomial, g_monomial)
                break
        if not new_monomial:
            # No monomial of matching degree was found
            h.add_monomial(monomial)
        else:
            h.add_monomial(new_monomial)

    # Add remaining unique monomials
    for monomial in g_monomials:
        h.add_monomial(monomial)

    return reduce_polynomial(h)

def multiply_polynomials(f, g):
    check_ordering([f, g])

    monomials = []
    for f_monomial in f.monomials:
        for g_monomial in g.monomials:
            monomials.append(multiply_monomials(f_monomial, g_monomial))

    return reduce_polynomial(Polynomial(monomials))

def zero_polynomial(f):
    for monomial in f.monomials:
        if monomial.coefficient != 0:
            return False

    return True

def reduce_polynomial(f):
    new_monomials = []
    for monomial in f.monomials:
        if monomial.coefficient != 0:
            new_monomials.append(monomial)

    return Polynomial(new_monomials)

def check_ordering(fs, order="lex"):
    incompatible_fs = filter(lambda x: x.order != "lex", fs)

    if incompatible_fs:
        raise Exception("Incompatible monomial orderings!")

if __name__ == "__main__":
    ## Eg. #1, Page 59
    ## xy^2 + 1
    #f = Polynomial([(Fraction(1), (1, 2)), (Fraction(1), (0, 0))])
    ## xy + 1
    #g = Polynomial([(Fraction(1), (1, 1)), (Fraction(1), (0, 0))])
    ## y + 1
    #h = Polynomial([(Fraction(1), (0, 1)), (Fraction(1), (0, 0))])

    # Eg. #2, Page 60
    # x^2y + xy^2 + y^2
    f = Polynomial([(Fraction(1), (2, 1)), (Fraction(1), (1, 2)), (Fraction(1), (0, 2))])
    # xy - 1
    g = Polynomial([(Fraction(1), (1, 1)), (Fraction(-1), (0, 0))])
    # y^2 - 1
    h = Polynomial([(Fraction(1), (0, 2)), (Fraction(-1), (0, 0))])

    result = divide_polynomials([g, h], f)
    print result
